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^P , Abstract 

G^ 

0^ ■ It is generally believed that the dynamics of simple fluids can be 

considered to be chaotic, at least to the extent that they can be mod- 
^-^' eled as classical systems of particles interacting with short range, re- 

' O I pulsive forces. Here we give a brief introduction to those parts of 

Q ■ chaos theory that are relevant for understanding some features of non- 

C^ I equilibrium processes in fluids. We introduce the notions of Lya- 

i-G ■ punov exponents, Kolmogorov-Sinai entropy and related quantities us- 

ing some simple low-dimensional systems as "toy" models of the more 
complicated systems encountered in the study of fluids. We then show 
how familiar methods used in the kinetic theory of gases can be em- 
V^ ' ployed for explicit, analytical calculations of the largest Lyapunov ex- 

C^ I ponent and KS entropy for dilute gases composed of hard spheres in 

d dimensions. We conclude with a brief discussion of interesting, open 
problems. 

1 Introduction 

We consider here a classical many-particle system, a gas of hard spheres or 
of hard disks. Our principal concern will be to develop methods by means 
of which we can understand and calculate the properties of such gases as 
chaotic dynamical systems. It is, of course, well known that to describe the 
macroscopic, equilibrium properties of such gases, we can easily dispense 
with any knowledge of most of the dynamical properties of the particles of 
which the gas is composed. That is, one can use thermodynamics and equi- 
librium statistical mechanics, i.e. statistical thermodynamics, to describe 



the relevant equilibrium properties of the gas. All of the relevant micro- 
scopic properties of the system needed for statistical thermodynamics are 
contained in the partition sum, which is defined in terms of the Hamiltonian 
of the system. The partition function is based on a probability measure on 
phase space. The macroscopic properties are simply related to averages with 
respect to this measure, of certain microscopic expressions. Of course it is 
far from trivial to compute these averages for anything like a real physical 
system. 

Our interest here, though, is to consider a gas as a mechanical system 
and to understand its behavior in time, rather than its equilibrium proper- 
ties, and to try to make quantitative statements about the motion of the 
trajectories of the phase points that describe the gas, in the usual 2Nd- 
dimensional phase-space, F-space, of the system, where A^ is the number of 
particles, d the number of the spatial dimensions of the system, d = 2 or 3, 
and the phase-space has dN spatial coordinates and dN momentum coordi- 
nates. We take the particles to be identical hard spheres or hard disks, each 
of mass m and diameter a. When we wish to describe the typical or aver- 
age properties of the system, we must start with the specification of some 
useful probability measure, with respect to which averages can be defined. 
Any dynamical system, therefore, consists of: 1) a space F, 2) a measure 
/i(j4), A (ZT, and 3) a transformation S" : F — > F. We will see that the dy- 
namical viewpoint can explain some features of macroscopic systems from 
their microscopic behavior. The explanations can be followed most easily 
in dynamical systems of very low dimensionality. However, even in simple 
low dimensional systems, dynamics may become so complicated that is is 
effectively impossible to follow the dynamics for long time, starting from 
a typical initial point, and we will be forced to consider typical behaviors 
using some appropriate probability measure. Our interest will be focused on 
chaotic systems which have the property that any uncertainty in the spec- 
ification of the exact initial state of the system will grow exponentially in 
time, to the point where the future of a phase-space point can no longer be 
predicted to within a reasonable accuracy[|l[|. But we can still say something 
about probabilities. 

It turns out that there is a close connection between the chaoticity of 
the system and issues like irreversibility on the macroscopic level and, for 
a gas of particles that interact with short-range forces, the validity of ki- 
netic theory [0]. This connection will first be outlined in section |2|, for low 
dimensional systems. A more extensive treatment can be found in Ref. y]. 
In section ^ we return to a high dimensional system in the form of a hard 
sphere gas in equilibrium. At low densities we can use kinetic theory to 
calculate a measure of chaoticity called the largest Lyapunov exponent. In 
section |^ another chaotic characteristic of this system is calculated using 
kinetic theory: the Kolmogorov-Sinai entropy. In section ^ we make some 
concluding remarks and present some open questions. 



2 Dynamical Systems 

The standard approaches to the theory of non-equihbrium processes in flu- 
ids are based on three foundational pihars: (1) The identification of the 
macroscopic quantities of physical interest as averages of microscopic quan- 
tities over an appropriate ensemble of similarly prepared systems; (2) The 
use of the Liouville equation, either in its classical or in its quantum me- 
chanical version, to compute the time evolution of the ensemble distribu- 
tion function; and (3) The utilization of some kind of physically reasonable 
factorization assumption for the ensemble distribution function in order to 
transform the Liouville equation into a tractable equation whose solution 
can be used to make quantitative statements about the macroscopic quan- 
tities. Such a procedure is followed in the derivation of the Navier-Stokes 
fluid dynamic equations from the Liouville equation Q] for general fluids, and 
in the derivation of the Boltzmann transport equation, and its extensions 
to higher densities, from the Liouville equation for dilute and moderately 
dense gases. More phenomenological approaches to irreversible behavior in 
fluids often depend on explicit stochastic assumptions about the underlying 
dynamical processes taking place in the fluid |p. 

While they are of the highest importance for the development of theories 
of irreversible processes in fluids, both of these approaches to irreversible be- 
havior leave the answers to some fundamental questions obscure. In partic- 
ular, these approaches offer only qualitative insights into the reasons for the 
validity of the stochastic assumptions imbedded in these various procedures 
- either through factorization assumptions, which in essence, are statements 
about correlations and probabilities - or through the replacement of the 
exact dynamics by a stochastic, Langevin-type, dynamics. Further, while 
the approaches outlined above do predict an approach to an equilibrium 
state, under the proper physical conditions, and do provide experimentally 
veriflable statements about the approach to equilibrium, they do not give a 
complete picture of why the system approaches an equilibrium state, based 
upon the underlying microscopic dynamics. The general arguments for the 
use of stochastic methods are based upon the randomness of the microscopic 
motions of the particles but are not much more speciflc than that. The pic- 
ture that we have of the approach to equilibrium is generally based on the 
idea that the local averages of conserved mechanical quantities, such as mass, 
momentum, and energy, change very slowly in time compared to the local 
averages of nonconserved quantities. Thus the macroscopic behavior will be 
dominated by the slowest variables in the system, the local conserved quan- 
tities, and the equilibrium state will be achieved when these quantities have 
reached steady, homogeneous values. This picture, suggested by solutions 
of the Boltzmann equation, has led to important advances in the theory of 
fluids, among others, to mode-coupling theory. What is missing from it is a 
basic understanding of the necessary (or sufficient) properties of the inter- 



molecular potential for the system to approach an equilibrium state, as well 
as an understanding of the properties of the trajectories of the system, and 
the evolution of measures in F-space, that are responsible for the approach 
to equilibrium states and, perhaps, under different boundary conditions, to 
more complicated, but interesting non-equilibrium steady states. 

The application of ideas from dynamical systems theory to non-equilibrium 
statistical mechanics allows us to make some progress in resolving the issues 
described above. The application of ideas from chaos theory, in particular, 
enables us to make some quantitative statements about the type and degree 
of randomness of a dynamical system, even of large systems typically treated 
by statistical mechanics. It also allows us to describe equilibrium and non- 
equilibrium states of a system in terms of probability measures defined in 
F-space, and in terms of the time evolution of these measures. Moreover, 
there are interesting and unexpected connections between the macroscopic 
transport coefficients that describe the approach of a fluid to equilibrium, 
and microscopic quantities that describe the chaotic behavior of a fluid, con- 
sidered as a large, dynamical system. In this section we will outline some 
of these rather new ideas, and illustrate their applications to statistical me- 
chanics by seeing how they work for systems of low dimensions and then 
generalizing them, when possible, to higher dimensional systems. We be- 
gin with a very simple two-dimensional reversible system, the baker's map, 
which exhibits many of the features we would like to see in more general, 
higher dimensional systems. 

2.1 The baker's map 

The simplest example of a reversible system with chaotic dynamics is prob- 
ably the baker's map. Here we consider a two-dimensional phase space on 
a unit square. That is, F = (x,y);0 < x,y < 1. The map, B, operates only 
at discrete time steps, and moves points {x,y) to B(x,y) = {x',y') given by 



= ((ylm) f-V2<.<i. (1) 

This map is illustrated in Fig. Q. It is immediately clear that this map 
possesses an inverse, B~^, given by 

= [ Y^ J for < y < 1/2 and 





%y!^i'^ ) for 1/2 < y < 1. (2) 




((),()) 1 1- 1 2 X (0,0) 

Figure 1: The baker's transformation. 




The baker's map is clearly area-preserving, and it is time-reversible in 
the sense that the transformation T : {x,y) -^ (1 — y,l — x) serves as a 
time reversal transformation for this map, such that T o B o T = B~^, and 
T o T = 1, where 1 is the unit operator. 

Now we regard the unit square as a "toy" phase-space. The dynamics 
of the baker's map in this phase-space has the following properties: 1) Con- 
sider two infinitesimally separated points. Unless they have exactly the same 
x-coordinates, or the same y-coordinates, the images of these two points 
under several applications, called iterations, of the map B will cause the 
x-components to separate exponentially with the number of applications, 
with an exponent of A+ = In 2, whereas the y-components will converge 
exponentially to a common value with an exponent of A_ = — In 2. The ex- 
ponents, A-t, characterizing exponential separation or convergence of points 
in phase-space are called Lyapunov exponents. The directions in which the 
points converge exponentially, in this case just the y-direction, are called 
stable directions, and the directions in which they separate exponentially, in 
this case the x-direction, are called unstable directions. The fact that the 
Lyapunov exponents sum to zero is a simple consequence of the area pre- 
serving property of the baker's map, as can easily be seen by considering the 
evolution, with successive iterations, of the small rectangle with corners at 
{x, y), {x + 6x, y), {x, y + 5y), {x + 6x, y + Sy) under the baker map, B. This 
rectangle has constant area, but grows exponentially long in the x-direction, 
and exponentially thin in the y-direction. Sooner or later it gets stretched 
and folded in such a way that on a coarse grained scale, the unit square is 
covered uniformly. We will describe this behavior on a coarse-grained scale, 
by saying that the distribution of points becomes weakly uniform. It is im- 
portant to note that the projection of the small rectangle onto the x-axis 
will be uniform in a time, n„ on the order of 



n„ 



InSx 
"A~' 



(3) 



where A-|- = In 2, is the positive Lyapunov exponent for the baker's map. 
That is, the projection of the small rectangle on the unstable direction be- 
comes uniform much sooner than the distribution of points on the entire 
unit square becomes weakly uniform. 





2.2 The Arnold cat map and hyperbolic systems 

A different map with a similar dynamical behavior, i.e., area preserving, 
with exponentially separating and converging trajectories characterized by 
positive and negative Lyapunov exponents, is provided by the Arnold cat 
map, T{x,y), illustrated in Fig. |2[ and given by 

mod 1. (4) 

The area preserving property is guaranteed by the fact that the matrix rep- 
resentation of T has unit determinant, the integer coefficients of the matrix 
together with the mod 1 condition implies that the unit square, or more 
properly, the unit torus, is mapped smoothly onto itself by T. Such area 
preserving maps with integer coefficients are called toral automorphisms. 
The Arnold cat map also has stable and unstable directions associated with 
positive and negative Lyapunov exponents given by X± = ln[(3 ± ^5)/2]. 
In a way much like the baker's map, a small region of the unit square will 
be stretched and squeezed under the iterated action of the cat map, with 
projections on both the x and y-axes becoming uniform on a similar time 
scale as in Eq. (^), and with the distribution of points becoming weakly 
uniform on a longer time scale. 

The baker's map and the Arnold cat map are two simple examples of 
what are called hyperbolic dynamical systems. Briefly, and somewhat loosely 
stated, hyperbolic dynamical systems are defined by the action of some 
dynamical transformation, S on a phase-space F, such that: (a) one can 
identify stable and unstable directions in F, under the action of S, with 
negative and positive Lyapunov exponents, all bounded away from zero; 
(b) the stable and unstable manifolds (lines, surfaces, etc.) are continuous 
functions of the variables that define the phase space, and when the two 
manifolds intersect, they do so transversely; (c) the system is transitive, 
i.e., there exists some trajectory in the phase-space that is dense on the 
phase-space; and (d) for maps, i.e., for dynamical systems where S acts only 
at discrete times, there are no directions in phase-space with a Lyapunov 
exponent of zero, while for "flows", i.e. systems where S depends upon a 
continuous time parameter, the only direction in F with a zero Lyapunov 
exponent is the direction of the flow. Clearly the baker's map and the Arnold 
cat map are hyperbolic maps. A typical flow that one might examine for 
hyperbolicity is the motion of a phase point on surfaces of constant energy 
for a system of interacting particles. 

2.3 Ergodic and mixing systems 

The baker's map and the Arnold cat map are also examples of dynami- 
cal systems which are ergodic and mixing. Ergodicity is the property of a 
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Figure 2: The Arnold cat map. 



dynamical system that the time average of any integrable function of the 
phase-space variables will be equal to the ensemble average of this function, 
the average taken with respect to an appropriate time translation invariant, 
equilibrium measure. That is, if /(T) is an integrable function, then 



1 "~^ r 

iim-^/(s^r)=//(r)/.(dr), 

j=o 



(5) 



where fJ,{A) is a measure that is invariant, i.e., fJ-{A) = fi{S^^A) for any 
non-trivial set A, and ergodic, meaning that it is impossible to divide the 
whole phase-space into two invariant sets, each of positive measure]^. It is 
generally assumed, but not always proved, that our systems possess a unique 
ergodic measure. Students of statistical mechanics will naturally associate 
the idea of an ergodic system with the name of Boltzmann who used this 

^Eq. (|5|) doesn't have to hold for all points F, as long as the set of points violating it 
has measure zero, with respect to the measure in the definition of the dynamical system. 



idea to base equilibrium statistical mechanics on the laws of mechanicsR. 

Much of equilibrium statistical mechanics can be based on the laws of 
large numbers, and strict ergodicity, in the sense of Boltzmann, is not that 
essential. However, non-equilibrium statistical mechanics requires some deep 
underpinnings from mechanics, or from the theory of stochastic processes. 
Here we take the point of view that Hamiltonian mechanics is all that is 
needed, but that is certainly not the only possible point of view. For non- 
equilibrium statistical mechanics, it is useful to explore an idea of Gibbs, 
which is called the mixing property of a dynamical system. Mixing systems 
are always ergodic, but the reverse is not always true. To define a mixing 
system, we consider two arbitrary sets in the phase-space, A and B, say, 
both of nonzero measure, and the evolution of the set A in time. Suppose 
after n iterations of the map S the set A has moved to S"A, then the system 
is mixing if 

n^^ ll{B) /i(r) ' ^ ' 

where //(F) is the measure of the entire phase-space, such as the unit square 
for the baker's map or the cat map, or the constant energy surface for a 
more general system. The mixing condition simply means that the time 
evolution of a set in phase-space is such that, in a coarse grained sense, it 
gets uniformly distributed, with respect to the measure /i, over the entire 
phase-space. It can be proved rather easily that for a mixing system, non- 
equilibrium averages of integrable functions / will approach their equilibrium 
values in the course of time. 

2.4 The approach to equihbrium 

A nice illustration of the approach to equilibrium, as provided by baker 
or cat maps, is to consider the behavior in time of reduced distribution 
functions. That is, if we think of the unit square, again, as a phase space, 
then we can define a phase-space distribution function, /9n(x, y), as a function 
of the number of iterations of the map, n, and the coordinates x and y. 
The phase-space distribution function satisfies a discrete-time version of the 
Liouville equation, which is a form of the Probenius-Perron equation for area 
preserving maps. The appropriate equation is 

/9„(x,y) =/9„_i(B"^(x,y)), (7) 

which, written out in full detail, becomes 

Pn{x,y) = /0„-i(x/2,22/) for < y < 1/2 

= pn-i{{x + l)/2,2y-l) forl/2<y<l. (8) 



Traditionally, statistical mechanical systems were called ergodic if they are ergodic 
under Hamiltonian flow with the Liouville measure on the energy shell. 



A similar, but somewhat more complicated equation could be given for the 
Arnold cat map, but we will not use it here. 

For these simple two dimensional models, a reduced distribution func- 
tion is obtained by integrating the distribution function over one of the two 
phase-space variables, x or y. This integration is motivated by the fact that 
for a system of A^ particles, we are not particularly interested in the full 
A^-particle distribution function, but rather in the one or two-particle distri- 
bution functions that can be used to evaluate the macroscopic quantities of 
interest, such as mass, momentum, and energy densities. Since our simple 
maps have only two coordinates, we can only consider the very simple case 
where a reduced distribution function is obtained by integrating over one 
of the coordinates. For the baker's map we will construct the distribution 
function for the density of points in the x-direction, for reasons that will be- 
come clear as we proceed. That is, we define a reduced distribution function 
Wn{x) by 

Wn{x)= f dypn{x,y). (9) 

Jo 

Using Eq. (g), we can easily obtain a difference equation for Wn{x), as 



Wn{x) = \ 



W^n-l(f) + ^n-l(^; 



(10) 



This equation is, among other things, the Frobenius- Perron equation for 
the one-dimensional map, x' = 2x (mod 1), on the interval (0,1). What is 
more important here, though, is that except for very special initial values 
for Wo{x), such as Dirac delta functions on the periodic points of the map, 
Wn{x) approaches a constant, independent of x, as n ^ oo. This may be 
proved in a number of ways, but may be understood most simply by just 
drawing some possible functional forms for VFo(x) and follow what happens 
to them after a few iterations of Eq. ([lO| ) . A standard procedure is to make 
a Fourier expansion of Wo{x), and to notice that only the constant term 
remains as the number of iterations gets large. The approach to equilibrium 
in this simple system can be associated with the properties of the expanding 
manifold in our simple two-dimensional phase-space. Because of the stretch- 
ing of regions in phase-space in the unstable directions, functions defined on 
the unstable manifold will get "smoothed out" in the course of time, much 
the same way that a ball of dough gets smoother and smoother along the 
direction that the baker stretches it. The initial wrinkles in the phase-space 
distribution function will not get smoothed out along the stable direction, 
on the contrary, they typically will get more and more wrinkled as the sys- 
tem evolves. From these considerations we can see that the integration of 
the phase-space distribution over the stable direction in Eq. (^) was not 
chosen accidentally; had we integrated over x instead, we would not have 
obtained an equation with a nice equilibrium solution as n ^ oo. In fact, a 
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Figure 3: The projection, W(x,n), onto the x-coordinate, of the phase- 
space distribution function for the Arnold cat map. 

typical initial distribution will become smooth in the expanding directions 
but very striated in the contracting directions. However, eventually it will 
look uniform on a coarse grained scale, consistent with the mixing behavior 
of the baker's map. 

The connection between the approach to equilibrium and the expanding 
direction of a measure preserving map can be further explored by consider- 
ing the Arnold cat map. Here, the expanding direction is along a line that is 
not aligned along either of the coordinate axes. One would expect that for 
this model a projection of the phase-space distribution function along either 
the X-axis, or the y-axis, would approach an equilibrium value. That this is 
so can be seen from a simple computer calculation. We start with a phase- 
space distribution that is concentrated in a small region < x,y < 0.1. We 
then follow the evolution in time of x and y projections of the distribution 
function. In Figs. |3| and ^ we can easily see that these distribution approach 
constant values after three or four iterations of the map, much before the en- 
tire phase-space distribution is smooth, even on a reasonably coarse grained 
scale, which for this arrangement takes eight to ten iterations. This obser- 
vation may be generalized: reduced distribution functions on some lower 
dimensional projection of phase-space will always become smooth under the 
dynamics, unless the projected space is entirely spanned by stable directions 
(hence is some subset of the stable manifold) . 
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Figure 4: The projection, G(y, n), onto the y-coordinate, of the phase-space 
distribution function for the Arnold cat map. 

Not only does one see an approach to an equilibrium distribution for the 
projected distribution functions for these maps, one also sees that a suitably 
defined Boltzmann iJ-function decreases monotonically as the number of 
iterations increases. This is illustrated in Fig. ^ for the Arnold cat map, for 
both projected distribution functions, starting from the same initial state 
as described above. The figure shows the if-function for both projections, 
as calculated on a computer. For the baker's map, we can easily show the 
monotonic decrease in the H function analytically. To do this we define the 
i?-function by 

H{n)= f dx Wn{x) In Wn{x). (11) 

Jo 

If we now use the recursion relation for Wn{x), Eq. (p!oD, we find that 

H{n + l)= I dxWn+i{x) InWn+iix) 
Jo 



= H(n). 



Wn{^) + Wn{^^: 



In 



[Wn{^) + Wni^) 



Wn{^)lnW^{^) + W^{^^)lnWn{^^: 



(12) 
The inequality in Eq. (^) follows from the fact that f[{a + b)/2] < [/(a) + 
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Figure 5: The Boltzmann i7-functions H\Y{n) and Hcin), obtained by 
using W{x, n) and G{y, n) respectively. 

f{b)]/2 if f{x) = xlnx. That is, a chord connecting two points on the curve 
xlna; hes above the curve. Thus we see that the H function decreases with 
time until Wn{x) becomes a constant. 

We conclude this discussion with a few remarks. For the baker's map 
and the Arnold cat map, admittedly "toy" models, highly simplified ver- 
sions of N particle systems, with almost trivial phase-spaces, we have been 
able to derive irreversible equations and i/-theorems with a minimum of 
assumptions. We have associated the approach to equilibrium of projected 
distribution functions with the existence of unstable manifolds for the dy- 
namics in phase-space, and the fact that the projection is not orthogonal to 
the unstable directions. In a more general context, such as a correspond- 
ing, but not yet possible, dynamical derivation of the Boltzmann transport 
equation, we would expect that the approach to equilibrium, seen here for 
baker and Arnold cat maps, would correspond to the approach to a local 
equilibrium state in the fluid. In such a local equilibrium state, the system 
has equilibrium values for density, temperature, and local mean velocity over 
distances of a few mean free paths. Then much slower hydrodynamic pro- 
cesses with a different kind of dynamics govern the approach to an overall 
equilibrium state for the entire fluid.^ 

^We wish to insert one word of caution about this picture. It seems clear from an 
examination of diffusion in some non-chaotic systems, such as the famous wind-tree model, 
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2.5 The Kolmogorov-Sinai entropy 

We next turn to a brief discussion of an important quantity that character- 
izes both deterministic chaos of the type we have been studying, as weU as 
Markov, stochastic processes. This quantity is called the Kolmogorov-Sinai 
(KS) entropy. For a deterministic system it measures the rate at which in- 
formation is gained about the initial state of a system. That is, suppose 
that we know that the initial phase point of our system is in some small 
region of F-space of dimension e on a side, and that we cannot resolve the 
location in F-space to any better precision. Consider now the evolution of 
this small volume element in F-space. For a system with non-zero Lyapunov 
exponents, this small volume will get exponentially stretched along the ex- 
panding directions. After some time this stretching will make some sides 
exponentially longer than the initial value e, typically of length eexp(Ait), 
where Aj is one of the positive Lyapunov exponents. Since we can resolve 
points in F-space to within a distance of e in any direction, we can now 
determine which of many small regions, of dimension e on a side, our sys- 
tem is in at time t. Then by inferring where this region came from in the 
initial volume, we learn more about the initial location of the phase point. 
Although it requires some careful analysis to prove, it is not difficult to 
imagine that there is some direct relation between the positive Lyapunov 
exponents and the KS entropy, hxs- In fact for a closed, hyperbolic system, 
Pesin has proved that the relation is as direct as our discussion above would 
imply, namely, 

hKS = E ^- (13) 

Ai>0 

A deep and interesting fact is that at least one way to prove Pesin's result 
(which we will not do here) depends on the fact that hyperbolic systems can 
be mapped onto Markov stochastic processes. That is, for such systems, 
with baker and Arnold cat maps as simple examples, one can represent the 
dynamics to within any arbitrary degree of precision, as a Markov process. 
Such Markov processes have a measure of their own information entropy, a 
quantity which measures the degree of uncertainty in the next outcome of 
the stochastic process. One can show that the information entropy, suitably 
defined, of the stochastic representation of a hyperbolic dynamical system 
is equal to the KS entropy of the system. This is one of the deep results of 
dynamical systems theory, which provides a firm mathematical basis for the 
correspondence of hyperbolic dynamical systems with Markov processes. It 
expresses the fact that from the point of view of mathematical analysis, at 



that chaoticity is sufficient, but not always necessary for understanding the approach to 
equilibrium of systems of many particles. However, in such non chaotic systems there 
often is a non-dynamical source of randomness, as in the random locations of scatterers in 
the wind tree model. This non-dynamical source of randomness is not needed to explain 
the approach to equilibrium in chaotic systems. 
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least, there is no real difference between a Newtonian, hyperbolic dynamical 
system with a finite KS entropy, and a Markov stochastic process with equal 
values of the information entropy. From the point of view of physics there is 
a big difference, of course. The central idea of the physical approach we take 
is to show that Newtonian dynamical systems are sufficiently hyperbolic to 
behave as if they were Markov stochastic systems and in consequence all 
important properties of Markov systems apply to dynamical systems, too. 

An important first step is to show that some useful models of physical 
systems have positive Lyapunov exponents and a finite, positive KS entropy 
per particle. In the next sections we will show how methods of statistical 
mechanics can be used to calculate Lyapunov exponents and KS entropies of 
some simple many-particle systems — gases of hard spheres in d dimensions, 
where d = 2,3. Before doing so, however, we briefly turn our attention 
to an application of the ideas of this section to non-equilibrium statistical 
mechanics. 

2.6 The escape-rate formalism for transport coefficients 

We conclude this section with a discussion of a formal relation between 
the transport coefficients that characterize hydrodynamic processes in fiuid 
systems, and the chaotic properties, such as Lyapunov exponents and KS 
entropies that characterize the underlying dynamical behavior of the fluid. 
The relation of interest here is called the "escape-rate" formula for transport 
coefficients and is due to Gaspard, Nicolis, and co-workers[P, Q. It applies 
to those fluids which can be considered to be classical, transitive, hyperbolic 
dynamical systems. While we do not know with certainty if any models of 
fluid systems satisfy this hyperbolicity requirement, we can suppose, as a 
working hypothesis, that generic fluid systems are well described as tran- 
sitive hyperbolic systems, and then explore the consequences that result. 
This hypothesis has been assumed by almost all workers in this field, but 
its most satisfactory articulation was provided by Cohen and Gallavotti in 
their study of non-equilibrium fluctuations in thermostatted systems [^Q 

To illustrate the escape-rate formula, we consider only the case of particle 
diffusion in an array of flxed scatterers, and refer to the literature for more 
general cases 0, ^. We suppose that a collection of particles is moving in a 
region R in space, and that there is also a collection of flxed scatterers placed 
in R, as illustrated in Fig. @. We suppose that the size of R is characterized 
by a length L which is much larger than the typical mean free path of the 
particles moving in R. For simplicity we suppose that the moving particles 



*It is customary to use the term Anosov system to describe a transitive hyperbolic 
system without singularities, such as the Arnold cat map. The class of Anosov systems 
does not include baker maps or hard sphere systems, since discontinuities of the map or 
flow, present in these systems, are not allowed by the Anosov condition. For this reason 
we prefer to generalize the chaotic hypothesis to include transitive, hyperbolic systems. 
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Figure 6: A slab geometry for diffusion in a system of moving particles is 
an array of fixed scatterers, with absorbing boundaries. 

do not interact with each other but only with the scatterers, and that the 
scatterers do not "trap" the moving particles in regions of microscopic or 
macroscopic size. In order to provide a non-equilibrium situation which 
would exhibit the diffusive properties of this arrangement, we suppose that 
R is surrounded by absorbing boundaries such that if a particle crosses the 
boundary of R from the interior, it is absorbed and lost to the system. 

Under most circumstances [|, the classical, macroscopic dynamics of this 
process is described by the diffusion equation 



dn{f, t) 

at 



DV'^n{f,t), 



(14) 



where n{f, t) is the density of moving particles at time t at point r, and D 
is the coefficient of diffusion of the moving particles for this system. The 
absorbing boundary conditions require that n(r, t) = on the boundaries. 
While the exact solution of Eq. ( [T^ ) depends on the geometry of the system, 
one can easily see that for long times, the total number of particles in the 
system decays with time as 



Nit) 



R 



dfn{f,t) ~ N{t = 0)exp[-DAt/L'^], 



(15) 



where A is a factor of order unity that depends on the geometry of the 
region R, D is the diffusion coefficient, and L is the characteristic size of 
R. We see that this macroscopic process is characterized by an exponential 
escape of particles from R, with an escape-rate, 'jmac = DA/L'^. There is a 
corresponding microscopic description of the escape process based upon the 



I.e., the surface area to volume of R scales to zero as L — » oo. 
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properties of the trajectories of the particles moving in the array of scatter- 
ers. There is a set of initial points in the phase space for the moving particles 
which are associated with trajectories that never leave the system for either 
the forward or time reversed motion. This set of points is called a repeller 
and denoted by TZ, and it is invariant, in the sense that any time-translation 
of this set of points is the set 7^ itself. Further, the repeller usually forms 
a fractal set in phase space with zero Lebesgue measure. Simple examples 
of repellers can be found in most books on chaos theory [Q, ^. It is possible 
to consider the dynamical properties of the trajectories on the repeller, and 
to define the Lyapunov exponents, Xi{TZ), and the KS entropy, hxsiT^), of 
such trajectories ||9|, |l^, 11|. For transitive hyperbolic systems, one finds that 



the sum of the positive Lyapunov exponents on the repeller is not equal to 
the KS entropy as would be the case for closed systems according to Pesin's 
theorem, but that the two quantities differ by an amount equal to the rate 
at which the other trajectories escape from the system, which we denote by 
7mic. That is 

-fmic = Y. ^^(^) - f^Ksin). (16) 

Ai>0 

Now we make the reasonable conjecture that the microscopic and the macro- 
scopic escape-rates are equal, which leads to an expression for the diffusion 
coefficient D given by 



1-2 

D = lim — 

L^oo A 



J2 A^(7^) - hKsin) 

X,>0 



(17) 



where we have taken the large system limit in order to remove terms of 
higher order in 1/L resulting from deviations of the actual dynamics from the 
diffusion law. In the event that the limit on the right hand side of Eq. ( [l7| ) 
exists, one has an expression for a macroscopic transport coefficient in terms 
of microscopic dynamical quantities. The escape-rate formalism has been 



applied by Gaspard and Baras|12] to determine the diffusion coefficient of a 
particle moving in a dense array of hard disk scatterers, where the centers 
of the scatterers are placed at the vertices of a triangular lattice, and by van 
Beijeren, Dorfman, and Latz, to determine the KS entropy on the repeller 
of a dilute, random Lorentz gas with hard disk or hard sphere scatterers p^]. 

3 Largest Lyapunov exponent of a gas of hard 
spheres at low density 

3.1 The hard sphere gas 

Often the calculation of chaotic characteristics of a system can only be done 
numerically. It would be preferable if one could at least find approximate 
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values for such quantities using more analytical methods, and thus gain 
some insight into the relevant physical processes. For the Lorentz Gas at 
low densities, this was done by Dorfman, Van Beijeren and others |15, ^, 



16, 17]. Here we present a calculation of the largest Lyapunov exponent for 
a system that is closer to a real gas, namely a dilute hard sphere gas. A 
brief presentation of this calculation can be found in Ref. flq] . 

We take A^ hard spheres in a volume V, in d dimensions. The diameter 
of the hard spheres is a, the reduced density is defined as the dimensionless 
number h = Na'^/V. We work in the thermodynamic limit N,V ^>- oo, 
keeping n fixed (but small). In this limit we need not be concerned with 
boundary conditions, but one may think of periodic boundary conditions, 
which have been used in Molecular Dynamics simulations to which we will 
eventually compare our results. 

The phase-space of the hard sphere gas consists of the positions {fj} and 
velocities {vi} of all N particles. To calculate the largest Lyapunov expo- 
nent, denoted here by A+, we need to consider two infinitesimally close 
trajectories in phase-space, F = {fi,vi, . . . ,fN,VN) and F + (5F = F + 
{6fi,Svi, . . . ,SfN,Svis[). The dynamics of the 6fi and 6vi is found from 
linearizing the dynamics of ri and Vi, which consists of a sequence of free 
flights and binary collisions. In free flight, there are continuous changes. 



n = 


= Vi 


Vi = 


= 


Sfi = 


= Svi 


5vi = 


= 0, 



(18) 

and at collisions the values of the two particles i and j change discontinuously 
according to: 

fj = rj 

Vj = Vj + {vij ■ a)a 

6fj = 6fj + {6fij ■ a)a 

6vj = 6vj + {6vij ■ a)a + Q{6fi — Sfj). (19) 

Primes are used to denote values right after the collision. Vij = Vi — Vj is 
the relative velocity and a = {fi— Tj)/a is the collision parameter. Q is the 
matrix [p!q] 

Q ^ [(^ • Vij)l + ^Vij] ■ [(^ • Vij)l - VijO-] ^20) 

a{a ■ Vij) 

The non-dotted products of vectors are dyadic products and 1 is the identity 
matrix. 

Our approach will be based on kinetic theory. We are concerned with 
the distribution of (r, v, 5r, 5v) as a function of time. For low densities, the 
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evolution of the distribution function / is described by a kinetic equation p9|]. 
This equation is based on the assumption that two colhding particles are 
uncorrelated, so that the probability of a collision between a particle with 
{fi,vi, 6fi, 6vi) and one with {f2,V2, 6f2,5v2) is proportional to the product 
f{fi,vi,6fi,6vi) f{f2,V2,Sr2,6v2). 

The kinetic equation can, unlike the ordinary Boltzmann equation, be 
expanded in powers of l/|lnn| to get the low density behavior of /, and 
thus of A+. We will take a different but roughly equivalent approach, we 
will derive the effective dynamics of the 6fi and 6vi for low densities, and 
use that to write down a kinetic equation. 

3.2 Low density dynamics — Clock model 

The main characteristics of the low density region is the typically long free 
flight time r of an individual particle between collisions, compared to the 
time it would take two transparent hard spheres to cross each other. If vq is 
the typical thermal velocity, the latter is (t/vq, while r ~ l/{voa N/V) = 
a/{vQh). Thus h is the small parameter. 

Just before a collision at time t the 5fi{t) will be 

5fi{t) = 5ri{to) + 5vi{to)Ti , (21) 

if to is the time of the previous collision and tj is the (large) free flight time 
t — to- Suppose that initially 6fi{to)/a and Svi{to)/vo are of the same order, 
then just before collision, 

6fi = Ti [5vi + 0{h)] . 

We insert this into the collision rules and neglect the terms of relative order 
0(n): 

5fj « TjSvj + {{jiSvi — Tj6vj) ■ a}cr 

6v- « Q{Ti5vi-Tj5vj). (22) 

Using Eq. (pO|), we see that Sf^/a and SvI/vq are both of order {6vi — 
6vj)/vQh, and are of the order of the ratio of the mean free time to the 
time it takes to traverse a particle diameter. For a dilute gas, this ratio is 
large and terms of relative order h can be neglected. In this way we have 
eliminated the 5fi from the 6vi dynamics. 

The neglected terms were of relative order n, relative to either Ti6vi or 
Tj6vj. These two are not necesarily of the same order. Now, if one of them 
is one or more orders of n higher than the other, we should also neglect it. 
If they are both of the same order, we should keep both. To know which 

®Of course one also has to demand that the particles should be a distance a apart. 
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terms to neglect, we have to keep track of the orders of n in 5vi. For that 
purpose, we define 

6vi=vo{ny'''ei, (23) 

where ||ej|| = 1. The number ki counts the number of orders of h and we 
will call this the clock value of particle i. The clock values are real numbers 
at this point, but later will be approximated by integers. Inserting Eq. ( |23| ) 
into Eq. (^), we can get the clock values k'^ and k'- after collision. Since, to 
leading order in density, 6vl and 6vj differ only in sign, k'^ = k', = k', with 



1 

Inn 



k' = — - — - In \\Q{Ti5vi — Tj5vj] 



Both Tj and tj are typically of order a /{v^h). This means that if ki > kj, 
we should neglect the term with 5vj, and if ki < kj, we should neglect the 
other term. This yields: 

k' = kD-\ r~^ln||QTDe£)||, 

— mn 

where D = i if ki > kj, and D = j otherwise. Particle D is called the 
dominant particle. Using the property td = 0{a /{vqu)) and the explicit 



form of Q from Eq. (20), one gets 

k' = kD + l + 0{-^). 

mn 

So far we have ignored the possibility that ki = kj. But the resulting 
correction in fact only contributes to the 0(1/ Inn) part. 

Differences in the number of collisions suffered by different particles, 
cause the clock values to not be all the same, even if they were so initially. 
They are indispensible for determining the magnitudes of postcollisional 
velocity deviations. But in fact no more is needed! For if we know how 
fast the clock values grow, we know the linear growth of In \\5v\\, which is 
precisely the largest Lyapunov exponent. So we define the clock speed as 

< k{t) > 
w = lim z ) 

in which < k > is the average clock value and v' is the average collision 
frequency.^] Because we extracted i>, which is 0{n), this clock speed is of 
order 1. 

The Lyapunov exponent is related to w via 

A_i_ = —wulnh. 



'In a hard sphere gas in d dimensions, u — ^!^ f \/ T^^ ^- 
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The clock speed w will be calculated in an expansion in 1/| In h\ . The leading 
order for low density, calculated in the next section, is a non-trivial constant. 
This behavior of A_|_ as a function of density had been conjectured already 
by Krylov pOJl , however without a nontrivial prefactor w. A first estimate of 
w was made by Stoddard and Ford pl[| , who got w = \nN. This would mean 
that there is no thermodynamic limit for the Lyapunov exponent. Some nu- 
merical simulations |2^] have been interpreted as supporting the logarithmic 
divergence with A^, but with a prefactor much smaller than 1. However, we 
will find a finite thermodynamic limit for w and show that, even in a mean- 
field approach that fully ignores local density fluctuations, it approaches this 
thermodynamic value so slowly that in the range of particle numbers acces- 
sible to simulations one could not distiguish between saturation or slow but 
steady increase. 

3.3 Kinetic approach 

For low densities, we may describe the effective dynamics of the clock values 

by 

k'- = kj = max{ki, kj) + l, (24) 

where the collision pairs {i,j) are chosen completely randomly with Poisson 
distributed collision times. The model with this dynamics we will call the 
clock model. For simplicity we will consider integer clock values only, though 
this restriction is by no means necessary or important. The clock speed w 
found in this model gives the leading term in the density expansion of the 
Lyapunov exponent: 

X+ = wiy[-lnn + 0{l)]. (25) 

A distribution function /^(t) will denote the fraction of particles having 
clock value k at time t. From the dynamics specified above we can derive an 
equation for the distribution function fk{t) of clock values. We expect the 
clock values to grow linearly with time. If they all grow at the same rate, 
we have 

fk{t) = g{k-w9t), 

with w as defined before, because 

-. oo -. oo 

lim ^ y^ g{k — wi>t)k = lim ^ y^ g(x)(x + wvt) = w. 

t^oo ut , ^-^ t^oo i/t ^-^ 

k = — 00 X = -rX) 

So once we have a kinetic equation for fkit), we will look for these propa- 
gating solutions. 

Consider the contributions to -^ from collisions in which a clock value 
k is lost. In each collision where a particle with k enters, the k gets lost, 
so the fraction of particles with k decreases at a rate yfk{t) due to these 
processes. There are also processes which increase the fraction of particles 
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having k. For these the larger incoming clock value should be A: — 1. We 
have to distiguish collisions with equal incoming clock values, both k — 1, 
from ones with different incoming clock values k — 1 and I. The rate for the 
latter type of collisions is i>fk-i Yl,i=-oo /'• -'^°'^ ^^^ former ones the rate is 
only 1/2 ^'/|_x, because the two particles are drawn from the same fraction 
/fc_i, and it doesn't matter in which order they are picked. In either case 
the number of particles with clock value k increases by two, so we get the 
kinetic equation: 

^^'^'^--A(t) + /|_, + 2A_i Y. /^ 



dut 



1= 



We scale out v by defining a new time variable r = vt. The equation is 
simplified further by replacing the fraction /^ of particles having clock value 
/c, by the fraction Ck of particles that have clock values k or less: 

Ck{r) = Y. M^)- 

l = — OD 

The kinetic equation then takes the short form: 

-Ck + Cl_,. (26) 



^C'fc ^ , ^2 



dr 

3.4 Front propagation 

Let us investigate the solutions of Eq. (|26|). We see that, given Ck-i and 
the initial value of C^, we can solve for C^: 



Ckir) 



Jo 



If Cfc = initially, it will remain zero, because < Ck-i < Ck- We take 
initial conditions such that Cfc<i(0) = and Cfc>o = 1, i.e. all particles have 
k = \. The solutions Ck are all polynomials in exp(— r), but with increasing 
k the order grows exponentially. Nonetheless, we calculated the Ck up to 
k = 2)2 using a computer program to handle the analytic manipulations. 
The Cfc's are plotted in Fig. | for r = 0, 2, 4, 6 and 10. 

We see in Fig. ^ that the initial distribution changes to some smoother 
shape, and moves to the right. After a while the shape seems to stay con- 



stant. We can now view Eq. (26) as describing a propagating front: C = 
is a stable phase and C = 1 is an unstable phase. On the left we have the 
stable phase, on the right, the unstable phase and in between is the interme- 
diate region called the front , that propagates to the right, into the unstable 
phase. The velocity at which it moves to the right is w. From Fig. ^ we see 
that for r = 10, the speed is still increasing, and is about 3.8. 
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Figure 7: Time evolution of the clock value distribution, 
values are 0. 



All initial clock 



Front propagation into an unstable phase comes in two flavors |23|| . In 
general, the instability of the unstable phase sets a velocity w* by which 
small perturbations propagate to the right. This velocity w* is determined 
from the linearized equation describing the front propagation around the 
unstable phase. For pulled fronts w* is the asymptotic velocity of any solu- 
tion with an initial shape that is sufficiently steep. If, however, no solution 
of the full non-linear front equation with this velocity exists, or if it is un- 
stable with respect to some nonlinear perturbation, the velocity is set by 
these non-linearities and the front is called pushed. In that case the velocity 
is higher than w*. We will assume that in our case the front is pulled. As 
we do not know of a general criterion for deciding whether a front is pushed 
or pulled, we will use the results of computer simulations for the validation 
of our assumption. 

We want to find a solution to Eq. (p6|) of the form 



CkiT)=F{k-WT)=F{x) 
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which means 



w- 



dF{x) 
dx 



-F{x)+F'^{x- 1) 



(27) 



Now C is a positive, increasing function of k, bounded between and 1. 
The function F should have all these properties too. So we are looking for 
a solution of a form such as depicted in Fig. |8[ 

1 




Figure 8: Clock value distribution as a propagating front 

For a pulled front, we have to investigate the Eq. (^) linearized around 
the unstable phase. Writing F = 1 — A, we get 

dA{x) 



w- 



dx 



-A{x) + 2A{x-l) + 0{A' 



(28) 



Given w, the asymptotic solution of Eq. (^) is given by a sum of exponentials 1 24] 

A(x)=^^ie-^'^ (29) 

i 

The possible values of 7j are found by inserting exp(— 72;) into the linearized 
equation (in case of degeneracies, Ai should be replaced by a polynomial in 
x). This gives w; as a function of 7: 



w{j) = {2e^ - l)/7. 



(30) 



The 7j may be complex, and, for given w, there are infinitely many of them. 
However, we know that A should be monotonic, so the most slowly decaying 
term in the sum should have a real positive 7. 

From the plot of Eq. ( pO|) in Fig. P, we see that there are three cases: 
If w is larger than some critical w* , there are two 7's, which are real and 
positive. If w = w* these two become degenerate and for w < w* they 
become complex. But for the slowest term to have complex 7 was not 
allowed by the condition of monotonicity of F: asymptotically the function 
would oscillate. So we conclude that, for the the Ansatz of a propagating 
front to work, the velocity should be at least w*, the minimal w from Eq. (|30| ) 
for positive real 7. This value can be expressed in terms of Lambert's W 
function:^ 

w* = __.~\. « 4.31107... 



m^ 



Defined as W{x) exp[W{x)] = x, where the branch analytic in is meant. 
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Figure 9: Velocity versus asymptotic decay rate 7. 

and it is the velocity set by the instability that we mentioned before. So it 
is the velocity w we were after. Its value is compatible with the estimate 
3.8 from Fig. |^ (a lower value than 3.8 would not), which supports the 
assumption of a pulled front. This concludes the calculation of the leading 
order in Eq. (^) of the largest Lyapunov exponent in the infinite system 
limit within the framework of our clock model, but we still have to consider 
the effects of a finite number of particles. 



3.5 Large finite A^ effects 

The kinetic equation works well in the thermodynamic limit, but to compare 
our results with simulations we need to correct for the effects of the finiteness 
of the number of particles. The clock model is very suitable for simple 
simulations. We take a set of A^ integer numbers {ki}. In each time step 
two of them are picked at random, and the collision rule in Eq. (|2^ ) is applied 
to them. We do this T times. The average number of collisions per particle 
is then 2T/N, because in each of the T collisions two particles are involved. 
An estimate for the clock speed is the average clock value J2i ^i/N divided 
by the average number of collisions per particle. For large T this approaches 
the clock speed W]\j, so 



1 



N 



WN 



lim — > kj 



2T 



This simulation is not even a bad chacterization of what happens with the 
clock values in the hard sphere gas, because at low densities there is very 
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little correlation between pairs of particles involved in subsequent collisions. 
It is very similar to the variant of the Direct Simulation Monte Carlo method 
used by Dellago and Posch|25] for the calculation of Lyapunov exponents in 



a "spatially homogeneous system". 

Clock speeds for different numbers of particles are plotted in Fig. ^j\ . 
Each point is determined from a single run with 2000 collisions per particle. 
The sums of all clock values at 20 times were fitted to a linear function. 
The slope gave wn, the error in the slope gave the error in wn- The errors 
are always less than 0.5 percent. One sees that wn increases slowly with 
increasing A^. It is possible that it saturates at the value of 4.311, but this 
cannot yet be seen even for half a million particles. 

A natural thought is that the N dependence of w is due to correlations 
between subsequent collisions: if two particles collide that already collided 
just before, or that had their clock values reset shortly before by parti- 
cles that were roughly synchronized already, the gain in clock value of the 
"slower" particle will be less than average. Hence the average clock speed 
is reduced. However, as we will see, the reduction in clock speed can be 
explained entirely on the basis of the linearized equation plus some simple 
bounds on its region of validity and this explanation does not seem to require 
any effects of correlated collisions. 

In the head of the distribution, the finiteness of the number of parti- 
cles becomes important when /^ = 0{1/N). For this reason, Brunet and 



Derrida|26| treat the finite A^ effects by the introduction of a cutoff e = 1/N 
in the equation, i.e. they modify the equation as soon as fk < £■ They 
distinguish three regions: one where the non-linear behavior is dominant, 
one where the linear equation holds, and one where the cutoff is effective. 
By glueing the solutions in these regions together, one obtains a new e- 
dependent front velocity. 

The introduction of a cutoff in something that is supposed to be a distri- 
bution function, which is an average over realizations, seems hard to justify. 
A more satisfactory way to obtain it, is to shift all clockvalues in each re- 
alization such that the particle with the largest clock value has x = 0, and 
then average the distribution function (this was also suggested by Kessler 
et a/[|^] ) . Then the cutoff occurs naturally. 

For X < there is a region where we can use the linearized equation. We 
consider the leading term in Eq. (pS|): 

A ~ -^e-'^o^. 

N 

The prefactor c/N is obtained from the fact that A = 0{1/N) at x = 
(c is order 1). The linear regime ends when A becomes of order 1, so for 
X = — ln( — )/7o. This is illustrated in Fig. IC. 



In contrast with the infinite system, now we only have to demand mono- 
tonicity and positivity from the solution of the linearized equation in this in- 
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Figure 10: The cutoff at the head of the distribution and the Hnear regime. 

terval of width ln(— )/7o. This means that a small imaginary value is allowed 
for the 7 in Eq. ( p9| ) with the smallest real part. We denote 7 = 7_r + ^77. 
The consequential oscillations, 

A = ao e~"'^^ cos(7/x + (jy), 

should not cause sign changes in the function or its derivative. The derivative 
can be sign definite if at most half a wavelength fits in the interval, i.e. if at 
most 

II 



IR-^ 



ln{N/c) 

For the leading behavior c can be set to 1. Positivity of A itself also poses 
an additional bound on 7/, but this can be shown to be a higher order effect. 
We expand ^(7) around its minimum at 79: 

wM = 1^(70 + (J7) ~ wijo) -\ — w"6'j -\ — w"'Sj + h.o. (higher orders), 

2 6 



where w" 
part gives 



d2 



d'^i 



^ 2 (70) and w'" = ■j-t(7o)- w is still real, so the imaginary 

= 6ji ^w"6jR + ^w"'[36jI - 57?] + h.o. I , 

where 5^ = S^r + 16 jj. To first order this says that w"6'yji = ^w"'d^]: the 
shift in the real part of 7 is higher order compared to the imaginary part. 
The new velocity is written as w = wq + 6w , and 6w is obtained from the 
real part of equation Eq. ( pO|) expanded to first order: 



dw 



\w" [5^1 - 57?} + \w"' {5^1 - 5^r5^]] + h.o. 



-w"5'~i'j + h.o. 



' 2\n^ {N/c) 



+ h.o.. 



(31) 



which coincides with Brunet and Derrida's result (when c is set to 1). They 
however needed to consider how the linear region connects to the others. 
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Figure 11: Effect of finite N on the clock speed. Diamonds are the simula- 
tions; the solid line is the prediction Eq. (|3^; the dashed hue is the asymp- 
totic value 4.311.... The inset shows the linear behavior of Xj ^wq — wn 
(plotted with error bars) as a function of IniV. The solid line in both is the 
fit to the form in Eq. (|3^ 



while we only need that A is of order 1 at the border with the non-linear 
region. From Eq. (pO|), one can show that 'w"jq = {wq — 1), so we find the 
result 

^ " ^ (32) 



WN = Wo 



21n^(iV/c) 



In Fig. 11 the simulation results are plotted together with a fit. We cal- 
culate dw"^'"^ = 1/ ^Jwq — wn, with wn taken from the simulations, which 
should be a linear function of In A^ for large N: 



5w 



-1/2 



N- 



6 In TV. 



(33) 



Indeed this behavior is seen in the inset of Fig. 11. According to Eq. (|3^ 
a = -Inc and b = 7r-V2/(tt'o - 1) = 0.247.... The fit to the dataf] 
shown in the inset yields b = 0.23 it 0.01, consistent with the theoretical 
predictionf^, and a = —0.07 ± 0.07. The value of c corresponding to this a 

®As the prediction is for large N, only points for A'' > 256 are used in the fit. 
^''By allowing changes in a and b simultaneously, one finds an appreciably larger range 
of acceptable values for b than the error in b indicates. 
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is c ss 1.07, so it is of order 1 as expected. 

The values from our simulations have been compared in Ref. [^] to 
molecular dynamics simulations of hard spheres, from which w is found 
from a fit of A+ to the form Eq. (P5|). The values agreed very well with the 
simulations. 

3.6 Further refinements 

Comparing the results from high precision simulations of the clock model 
on the one hand and of hard disk systems with exactly the same number 
of particles and equal collision frequency on the other hand, one finds that 
the dimensionless clock speed w in the latter is significantly higher than in 
the former; for a system of 10, 000 particles the clock model gives a w oi 
4.05 ± 0.01, whereas the corresponding hard disk value is 4.47 it 0.02. The 
cause of this difference is the velocity dependence of the collision frequency, 
which is an increasing function of speed. As a result particles in the head 
of the distribution tend to have a higher speed than average (a higher col- 
lision frequency enhances the clock value), resulting in a clock speed that 
is systematically higher than it would be in the case of a velocity indepen- 
dent collision frequency. The clock model does have a velocity independent 
collision frequency indeed, which explains the difference in w between this 
model and the hard disk system. 

An explicit calculation accounting for these effects and producing im- 
proved estimates for w in hard disk systems will appear soon[p^. 

4 Kolmogorov- Sinai entropy of a gas of hard spheres 
at low density 

In the previous section we reviewed the calculation by Van Zon and Van 
Beijeren of the largest Lyapunov exponent for a gas of hard disks at low 
densities. Here we will outline a related, but somewhat simpler calculation of 
the KS entropy of a gas of hard disks or of hard spheres at low density. This 
calculation, while not rigorous in a mathematical sense, is a strong indication 
of the chaotic behavior of such gases, and indicates that the chaoticity of a 
dilute hard disk or hard sphere gas persists in the thermodynamic limit. 

The starting point is the same as in the previous section. It is important 
to note that the dynamics of a hard sphere system is exactly described as 
free motion of the particles, punctuated by instantaneous, binary collisions 
between some pair of particles. To describe this motion and the quantities we 
need for the KS-entropy, we consider the dynamical behavior of the positions 
and coordinates of all the particles in the gas, {fi,vi,f2,V2, ■ ■ ■ ^tn^vn), as 
well as a set of deviation vectors which describe the motion of a pencil of 
nearby trajectories in phase-space, {5fi, 6vi, . . . , (5r/v, 6v]\f)- The equations 
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of motion for these quantities between the binary colhsions are given by 



Eq. (18) and the changes of these quantities at a colhsion between particles 



i and j are given in Eqs. (19,20) 



To proceed with the calculation of the KS-entropy, we will suppose that 
the Pesin formula holds, i.e., that the KS-entropy is the sum of the positive 
Lyapunov exponents. This sum can be obtained by considering the growth 
in time of the volume of the dA^-dimensional projection of an arbitrary, 
infinitesimal volume in the full 2(iA^-dimensional phase space. This observa- 
tion requires some explanation, as follows. The typical rate of separation of 
two arbitrary, but infinitesimally close, points in phase-space, for a hyper- 
bolic system, will be exponential and the rate will be given by the largest 
Lyapunov exponent. If we consider a typical two-dimensional, infinitesimal 
area in phase-space, then this area will grow exponentially with a rate deter- 
mined by the sum of the two largest Lyapunov exponents. In other words, 
the exponential growth of a typical infinitesimal n-dimensional subvolume 
in phase space is determined by the sum of the n largest Lyapunov expo- 
nents. Further, for a Hamiltonian system, the Lyapunov exponents come in 
plus-minus conjugate pairs, so that the sum of a conjugate pair of exponents 
is always zero. Consequently, the growth of an infinitesimal dA^-dimensional 
subvolume is determined by the sum of the dN non-negative Lyapunov expo- 
nents, and the volume of an infinitesimal 2(iA^-dimensional volume remains 
constant, in accord with Liouville's theorem. 

For the typical dN dimensional subvolume, we consider a volume formed 
by the projection of 2dN infinitesimal displacement vectors on velocity 
space, {6vi,6v2, ■ ■ ■ ,Sviy)- Given initial values for each of these vectors, 
as well as for all of the positions, r^, velocities, vi, and position deviation 
vectors, Sfi, we can follow their evolution in time and, in principle at least, 
determine the time dependence of an infinitesmal volume element in velocity 
space, which we denote as 5V{t). Then 

hxs = y] ^i = lim 7 In ' 



Ai>0 



i^oo t 6V{0) 

,. 1 /•* dlnSVir) 
lim - / dr- 



t^oo t Jo dr 

dt I 

The last line of Eq. (p^) is based on the assumption (still unproven) that 
a gas of hard spheres is ergodic, so that time averages can be replaced 
by equilibrium averages taken with respect to a microcanonical ensemble. 
Here this average is denoted by angular brackets. Since we are considering 
a volume element in velocity space, we can use the fact that the velocity 
displacement vectors do not change during the free flight motion of the 
particles between the collisions, but do change at a collision. Under these 
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circumstances, one can use elementary kinetic theory considerations to show 
that the final term on the right hand side of Eq. (^) is 

'±^) . (i:t„i.,.v\ ^!M_R^rr,,,svy (35) 



Here T12 is a binary collision operator, discussed in some detail in Refs. |2£ 



30 1, given by 

^12 = a'^-^ f da\vi2 ■ a\6in2 - aa)[Vail,2) - 1]. (36) 

Jvi2a<0 

In Eq. (|^), d is the number of spatial dimensions of the system, a again is 
the diameter of the spheres, fi2 = n — f2;vi2 = vi — V2, and the operator 
P5-(l,2) is a substitution operator which replaces the precollision values, 
ri,vi,f2,V2,Sfi,6vi,Sf2,6v2, by their post collision values, denoted with 



primes, given by Eqs. ( |19| , ^C| ). The unit vector a is an impact parameter, 
running in the direction of the line connecting the centers at collision and 
is integrated over a hemisphere corresponding to all allowed directions. 
At this point it is useful to express the precollision quantities 5rj as 

Sfi = Sfi{0) + TiSvi, (37) 

where (5fj(0) is the position displacement of particle i just after its previous 
collision, and tj is the time between the previous collision of particle i with 
some other particle, and the next collision involving particle i. To further 
simplify the expression for the KS-entropy, we now neglect, as in section ^, 
the initial displacement vectors, Sfi{0), when we calculate the change of the 
infinitesimal volume in velocity space at the (1,2) collision in Eq. (|35|). This 
turns out to be a serious approximation. It leads to the correct value for the 
leading density term in h^Sy at low density, but the first order correction 
to this term is obtained incorrectly in this approximation. This can be 
repaired, but at the cost of a much longer and intricate calculation which 
we will present elsewhere. 

If we insert 6fi = tiSvi and 6f2 = T2dv2 into the expression, Eq. (|l9|), 
for the post collision velocity deviations for particles 1 and 2, we find that 



6V' 
[Va{l, 2) - 1] In (5V = In — - = In | det M12I, (38) 

oV 



where 



M12 = l-2aa 



a 



(ui2 ■ a)l — vua + avi2 — ,^ ^"^ ^, da 

(wi2 • 0-) 



(39) 



In Eq. (|39D, T12 = {t\ + T2)/2 and 1 is the unit matrix. The determinants 
are easily evaluated. For d = 2, one finds 

|detM,2| = l + ^^i^, (40) 

(TCOS0 
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where is the angle of incidence in the 1,2 cohision and ranges over the 
values — 7r/2 < (j) < tt/2. A similar calculation for d = 3 shows that 

|detMi2| = l + ^^^(cos2</> + l) + f^^^^V- (41) 

a cos (p \ a / 

To obtain the leading term in the density, for low density gases, we keep the 
highest power of the time in each of the expressions for the determinant. 
Further, at low densities we can compute the ensemble averages appearing 



in Eq. (34) by ignoring possible pre-collision correlations between particles 
1 and 2, and using equilibrium values for the single-particle distribution 
functions appearing in the ensemble averages. In this way we find for d = 2 

hxs/N = — / dvi / dv2 / dri / dT2 / da\vi2 ■ o-\ x 

2n J J J J Jvi2-a<o 

xFi{vun)Fi{v2,T2) InTu + • • • , (42) 

where the normalized equilibrium single particle distribution functions, Fi{vi,Ti) 
are given, in d dimensions, by 

FM,n) = n(^)^/2z.(^l)e-^™^1/2e-(^^')-^ (43) 

Here n is the number density of the gas, /? = {kBT)~^ , where T is the 
gas temperature, and ks is Boltzmann's constant, z^(tTj) is the equilibrium 
collision frequency for a particle with velocity Vi. For two dimensions, the 
evaluation of the integrals leads directly to 

hKs/N='^[-Hna^) + ---] (44) 

where u = [{2'K^''^na) / {(3m)^''^] is the average collision frequency at equi- 
librium for a two-dimensional gas of hard disks. The terms left out are of 
higher order in the density. 

For three-dimensional gases, a parallel calculation leads to 

hKs/N = iy[-ln{-iTna^) + ■■■], (45) 

where for a gas of hard spheres (d = 3) the average collision frequency i^ = 
[(47r^'^n(7^)/(/3m)-'^'^]. These results are in excellent agreement with the nu- 
merical simulations of Dellago and Poschp^]. The higher order terms take, 
as mentioned earlier, considerably more work, and are discussed elsewhere pl|. 



5 Conclusions and Outlook 

In the previous sections we have reviewed some of the ideas that motivate 
the interest in the chaotic foundations of non-equilibrium processes in fluids. 
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We have provided an elementary discussion of transitive, hyperbolic dynam- 
ical systems such as the baker and the Arnold cat maps to illustrate some 
of the central notions and dynamical quantities. We then turned to the ap- 
plications of kinetic theory to compute the largest Lyapunov exponent and 
the KS entropy for a dilute gas of hard disks or hard spheres. The explicit 
results obtained by these methods are in good agreement with the results of 
computer simulations, and, apart from the corrections due to the velocity 
dependence of the collision frequency referred to at the end of section ^, 
represent the present state of the art in the analytical calculation of chaotic 
quantities for dilute gases with short range, repulsive forces. There are, 
however, still many open problems which need solving. Here we mention a 
few of them: 

1. We have a theory for the leading density behavior of the largest Lya- 
punov exponent for a dilute gas of disks or spheres. We also have some 
understanding of the number dependence of this quantity when we are 
not quite in the thermodynamic limit. However we do not know much 
about the higher density corrections to this exponent, nor do we know 
anything about the rest of the Lyapunov spectrum, other than the KS 
entropy per particle (in the thermodynamic limit). The determination 
of the complete spectrum would be quite an accomplishment. 

2. Recent results of Rom-Kedar and Turaev[32| imply that systems with 
short range repulsive forces, other than hard disks or spheres, may not 
be totally hyperbolic. Instead their phase spaces may have elliptic is- 
lands where the motion is not chaotic. It would be interesting to know 
first of all if there are any experimental or theoretical consequences of 
the existence of these elliptic regions for non equilibrium processes in 
real fluid systems and also whether such elliptic islands will persist for 
arbitrarily large energies. 

3. One important application of the methods described here is to the de- 
termination of the chaotic properties of thermostatted, driven systems 
for which a non-equilibrium steady state is reached and maintained. 
The general properties of such systems are described in the books of 



Hoover [p^] and of Evans and Morriss|34|, and a clear mathematical 
description has recently been given by Ruelle |]35[| . Of special interest 
is the perturbation of the Lyapunov spectrum produced by the ther- 
mostatted driving field. For dilute, random Lorentz gases it has been 
possible to use kinetic theory methods to determine the spectrum when 



the field is small[16|. It would be worthwhile to extend these results 



to larger fields and to gases where all of the particles are moving. 

4. The escape-rate formalism described here has two drawbacks: (a) It 
is not at all easy to describe the fractal repeller that forms in the 
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phase-space of a system with many degrees of freedom, (b) Even in 
those cases where the sum of the positive Lyapunov exponents on 
the repeller can be calculated analytically, the KS entropy is not yet 
directly accessible to analytic methods. Instead one has to use the 
transport coefficients and the sum of the positive Lyapunov exponents 
to infer the KS entropy of trajectories on the repeller. It would be 
very helpful to have a better understanding of the properties of high- 
dimensional repellers and to have an independent analytical means to 
compute the KS entropy of trajectories on the repeller. 

5. One of the main goals of current research in this area is to obtain, 
if possible, some deeper understanding of the dynamical basis of the 
laws of irreversible thermodynamics. For two dimensional diffusive 
models based upon the baker map, it has been possible to show that 
the laws of irreversible thermodynamics result from a careful analysis 
of the fractal structures that appear in the relevant phase spaces of 
these models when the systems are in non-equilibrium steady states 1 36, 
37, ^]. The main physical idea is that entropy is produced by the 
irreversible loss of information when changes are taking place in a 
system on very fine scales, beyond experimental resolution. However, 
as has been emphasized by other authors |39|], these models may be 
too simple and/or the thermostats considered may be too special to 
allow for any general conclusions to be drawn. This area of research 
is active and many issues remain to be understood. 

6. Our description of fluid systems as composed of classical particles in- 
teracting through repulsive, short range forces is certainly incomplete. 
Typical fluid systems are better modeled by short range forces with 
both attractive and repulsive regions. We do not yet know what effects 
a more careful analysis of interparticle forces will have on our picture 
of the chaotic behavior of fluids. An even more serious problem is con- 
nected with our use of classical mechanics to describe systems which 
are quantum mechanical in nature. We have almost no understanding 
of how to correctly obtain a quantum version of the classical chaotic 
picture of fluids, or even know for sure if such a thing is possible. 
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